% This program is the test of Saturated Vapor Pressure
% clc
% clear

var_name = 'gh';
% var_name = 'gd';
% var_name = 'us';
% var_name = 'vs';
it       = 16;

gravity = 9.80616;

nc_file  = '..\run\hope_output.n90.nc';

lon = ncread(nc_file,'lon'   );
lat = ncread(nc_file,'lat'   );
if strcmp(var_name,'gd')
    gh  = squeeze(ncread(nc_file,'gh',[1 1 1 1 it],[Inf Inf Inf Inf 1]));
    ghs = squeeze(ncread(nc_file,'ghs',[1 1 1    1],[Inf Inf Inf      1]));
    var = gh - ghs;
else
    var = squeeze(ncread(nc_file,var_name,[1 1 1 1 it],[Inf Inf Inf Inf 1]));
end

case_id = ncreadatt(nc_file,'/','case_id');
dx      = ncreadatt(nc_file,'/','dx'      );

coef = 1;
if(case_id==5)
    lon(lon>180) = lon(lon>180) - 360;
    xs = -180;
    xe = 180;
    if(strcmp(var_name,'gh'))
        coef = 1 / gravity;
    end
elseif(case_id==8)
    lon(lon<-270) = 360 - lon(lon<-270);
    lon(lon>90) = lon(lon>90) - 360;
    xs = -270;
    xe = 90;
else
    lon(lon<0) = 360 + lon((lon<0));
    xs = 0;
    xe = 360;
end

lon1d = double(reshape(lon,[],1));
lat1d = double(reshape(lat,[],1));
var1d = double(reshape(var,[],1));

[lon_plt,lat_plt,var_plt] = extend_var_plot(lon1d,lat1d,var1d);

res = dx;
x   = xs:res:xe;
y   = -90:res:90;

[lon2d,lat2d] = meshgrid(x,y);

var_plot = griddata(lon_plt,lat_plt,var_plt,lon2d,lat2d,'linear');

figure
plt = pcolor(lon2d,lat2d,var_plot*coef);
% shading interp
set(plt,'EdgeColor','None')
% set(gca,'CLim',[ 7.8494e+04,1.0350e+05])
% set(gca,'CLim',[ 4.9e4,5.9e4])
% set(gca,'CLim',[ -7e-9,7e-9])
colormap(jet)

if(case_id==5)
    if(strcmp(var_name,'gh'))
        hold on
        LevelList = 5050:50:5950;
        plt=contour(lon2d,lat2d,var_plot*coef,'LevelList',LevelList,'LineColor','k');
        %     plt=contour(lon2d,lat2d,var_plot,'LevelList',LevelList,'LineColor','k','ShowText','on');
    elseif(strcmp(var_name,'zonal_wind'))
        hold on
        LevelList = -30:10:60;
        plt=contour(lon2d,lat2d,var_plot,'LevelList',LevelList,'LineColor','k');
    elseif(strcmp(var_name,'meridional_wind'))
        hold on
        LevelList = -30:10:30;
        plt=contour(lon2d,lat2d,var_plot,'LevelList',LevelList,'LineColor','k');
    end
end

function [lon,lat,var] = extend_var_plot(lon_in,lat_in,var_in)
n = size(lon_in,1);
lon = zeros(3*n,1);
lat = zeros(3*n,1);
var = zeros(3*n,1);

lon_neg = lon_in - 360;
lon_pos = lon_in + 360;
lat_neg = lat_in;
lat_pos = lat_in;
var_neg = var_in;
var_pos = var_in;

lon(1+0*n:1*n,1) = lon_neg;
lon(1+1*n:2*n,1) = lon_in;
lon(1+2*n:3*n,1) = lon_pos;

lat(1+0*n:1*n,1) = lat_neg;
lat(1+1*n:2*n,1) = lat_in;
lat(1+2*n:3*n,1) = lat_pos;

var(1+0*n:1*n,1) = var_neg;
var(1+1*n:2*n,1) = var_in;
var(1+2*n:3*n,1) = var_pos;
end
